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Based on first- and second-order perturbation theory, we present a numerical study of the tempo¬ 
ral build-up and decay of unsteady acoustic fields and acoustic streaming flows actuated by vibrating 
walls in the transverse cross-sectional plane of a long straight microchannel under adiabatic con¬ 
ditions and assuming temperature-independent material parameters. The unsteady streaming flow 
is obtained by averaging the time-dependent velocity field over one oscillation period, and as time 
increases, it is shown to converge towards the well-known steady time-averaged solution calculated 
in the frequency domain. Scaling analysis reveals that the acoustic resonance builds up much faster 
than the acoustic streaming, implying that the radiation force may dominate over the drag force 
from streaming even for small particles. However, our numerical time-dependent analysis indicates 
that pulsed actuation does not reduce streaming significantly due to its slow decay. Our analysis 
also shows that for an acoustic resonance with a quality factor Q, the amplitude of the oscillating 
second-order velocity component is Q times larger than the usual second-order steady time-averaged 
velocity component. Consequently, the well-known criterion Vi -C Cg for the validity of the pertur¬ 
bation expansion is replaced by the more restrictive criterion Vi <C c^fQ. Our numerical model is 
available in the supplemental material in the form of Comsol model files and Matlab scripts. 

PACS numbers: 43.25.Nm, 43.20.Ks, 43.25.-|-y 


I. INTRODUCTION 

Acoustophoresis has successfully been used in many 
applications to manipulate particles in the size range 
from about 0.5 mm down to about 2 pm [I|. However, for 
smaller particles, the focusing by the acoustic radiation 
force is hindered by the drag force from the suspend¬ 
ing liquid, which is set in motion by the generation of 
an acoustic streaming flow [1, Q . This limits the use of 
acoustophoresis to manipulate sub-micrometer particles, 
relevant for application within medical, environmental, 
and food sciences, and it underlines a need for better un¬ 
derstanding of acoustic streaming and ways to circum¬ 
vent this limitation. 

The phenomenon of acoustic streami^ was first de¬ 
scribed theoretically by Lord Rayleigh [J in 1884, and 
has later been revisited, among others, by Schlicting Q, 
Nyborg Hamilton 0,11], Rednikov and Sadhal and 
Muller et al. 0 , to extend the fundamental treatment 
of the governing equations and to solve the equations for 
various open and closed geometries. 

Numerical methods have been applied in many stud¬ 
ies to predict the streaming phenomena observed in var¬ 
ious experiments. Muller et al. 0| developed a nu¬ 
merical scheme to solve the acoustic streaming in the 
cross section of a long straight microchannel, which re¬ 
solved the viscous acoustic boundary layers and described 
the interplay between the acoustic scattering force and 
the streaming-induced drag force on suspended particles. 
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This scheme was later extended to take into account the 
thermoviscous effects arising from the dependence of the 
fluid viscosity on the oscillating temperature field m- 
Lei et al. [12, [I2| have developed a numerical scheme 
based on the effective slip-velocity eq uations, originally 
proposed by Nyborg in 1953 [l3, Ulli which avoid the 
resolution of the thin boundary layers but still enable 
qualitative predictions of the three-dimensional stream¬ 
ing flows observed in microchannels and flat microfluidic 
chambers. To obtain quantitative results from such mod¬ 
els that does not resolve the acoustic boundary layers, 
Hahn et al. [13 developed an effective model to deter¬ 
mine the loss associated with the viscous stresses inside 
the thermoacoustic boundary layers, and apply this loss 
as an artificial bulk absorption coefficient. This enables 
the calculation of correct acoustic amplitudes, without 
resolving the thin acoustic boundary layers. Acoustic 
streaming in the cross section of a straight PDMS mi¬ 
crochannel exited by surface acoustic waves was studied 
numerically by Nama et al. [13 , describing the influence 
of the acoustically soft PDMS wall on the particle fo- 
cusability, and examining the possibilities of having two 
tunable counter-propagating surface acoustic waves. 

All of the above mentioned studies consider steady 
acoustic streaming flows. This is reasonable as the 
streaming flow reaches steady state typically in a few 
milliseconds, much faster than other relevant experimen¬ 
tal timescales. Furthermore, this allows for analytical 
solutions for the streaming velocity field in some special 
cases, and it makes it much easier to obtain numerical so¬ 
lutions. However, an experimental study by Hoyos and 
Castro [13 indicates that a pulsed actuation, instead of 
steady, can reduce the drag force from the streaming flow 
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relative to the radiation force and thus allowing the latter 
also to dominate manipulation of sub-micrometer parti¬ 
cles. This might provide an alternative method to the 
one proposed by Antfolk et al. [l^ , which used an almost 
square channel with overlapping resonances to create a 
streaming flow that did not counteract the focusing of 
sub-micrometer particles. 

To theoretically study the effects of a pulsed ultra¬ 
sound actuation, we need to solve the temporal evolu¬ 
tion of the acoustic resonance and streaming, which is 
the topic of the present work. Numerical solutions of the 
time-domain acoustic equations were used by Wang and 
Dual ( 2 ^ to calculate the time-averaged radiation force 
on a cylinder and the steady streaming around a cylin¬ 
der, both in a steady oscillating acoustic field. However, 
they did not present an analysis of the unsteady build-up 
of the acoustic resonance and the streaming flow. 

In this paper, we derive the second-order perturbation 
expansion of the time-dependent governing equations for 
the acoustic fields and streaming velocity, and solve them 
numerically for a long straight channel with acoustically 
hard walls and a rectangular cross section. The anal¬ 
ysis and results are divided into two sections: (1) A 
study of the transient build-up of the acoustic resonance 
and streaming from a initially quiescent state towards 
a steady oscillating acoustic field and a steady stream¬ 
ing flow. (2) An analysis of the response of the acous¬ 
tic field and the streaming flow to pulsed actuation, and 
quantifying whether this can lead to better focusability 
of sub-micrometer particles. 

In previous studies, such as only the peri¬ 

odic state of the acoustic resonance and the steady time- 
averaged streaming velocity are solved. When solving 
the time-dependent equations, we obtain a transient so¬ 
lution, which may also be averaged over one oscillation 
period to obtain an unsteady time-averaged solution. 


II. BASIC ADIABATIC ACOUSTIC THEORY 

In this section we derive the governing equations for 
the first- and second-order perturbations to unsteady 
acoustic fields in a compressible Newtonian fluid. We 
only consider acoustic perturbation in fluids, and treat 
the surrounding solid material as ideal rigid walls. Our 
treatment is based on textbook adiabatic acoustics [2l| 
and our previous study Ref. m of the purely periodic 
state. 


A. Adiabatic thermodynamics 

We employ the adiabatic approximation, which as¬ 
sumes that the entropy is conserved for any small fluid 
volume [l^. Consequently, the thermodynamic state of 
the fluid is described by only one independent thermo¬ 
dynamic variable, which we choose to be the pressure p. 
See Table U for parameter values. The changes dp in the 


TABLE I. lAPWS parameter values for pure water at ambient 
temperature 25 °C a nd p ressure 0.1013 MPa. For references 
see Sec. II-B in Ref. [Ill] . 


Parameter 

Symbol 

Value 

Unit 

Acoustic properties: 




Mass density 

Po 

9.971 X 

lO’^ 

kgm"® 

Speed of sound 


1.497 X 

10® 

m s~® 

Compressibility 

Ks 

4.477 X 

10-10 

Pa"® 

Transport properties: 




Shear viscosity 

V 

8.900 X 

10“'‘ 

Pa s 

Bulk viscosity 


2.485 X 

10“® 

Pa s 


density p from the equilibrium state is given by 


dp = pKg dp, 

where the isentropic compressibility is defined as 


dp 

dp 


1 

pc- 


,2 * 


( 1 ) 

( 2 ) 


B. Governing equations 


Mass conservation implies that the rate of change d^p 
of the density in a test volume with surface normal vector 
n is given by the influx (direction —n) of the mass current 
density pv. In differential form by Gauss’s theorem it is 

dtP = ^ ■[ - P'v]- (3a) 

Substituting d^p and Vp using Eq. (IT|), and dividing by 
p, the continuity equation (I3al) becomes 

Kgd^p = —V • V — KgV ■ Vp. (3b) 


Similarly, momentum conservation implies that the rate 
of change d^{pv) of the momentum density in the same 
test volume is given by the stress forces cr acting on the 
surface (with normal n), and the influx (direction —n) of 
the momentum current density pvv. In differential form, 
neglecting body forces, this becomes 


dt{pv) = V • [x - p 1 - pvv ], 
where the viscous stress tensor is defined as 

. 2 


T = rj 


Vu -I- (Vt))^ 


P 




{V -v) 1. 


(4a) 

(4b) 


Here 1 is the unit tensor and the superscript ”T” indi¬ 
cates tensor transposition. Using the continuity equation 
(I3al) . the momentum equation (I4al) is rewritten into the 
well-known Navier-Stokes form. 


pd^v = V • [t — p 1 ] — p{v ■ V)u, (4c) 


which is useful when solving problems in the time do¬ 
main. The equations (l3b)) and (l4cl) constitutes the non¬ 
linear governing equations which we will study by ap¬ 
plying the usual perturbation approach of small acoustic 
amplitudes. 
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C. First-order time-domain equations 


E. Periodic frequency-domain equations 


The homogeneous, isotropic, quiescent thermodynamic 
equilibrium state is taken to be the zeroth-order state in 
the acoustic perturbation expansion. Following standard 
perturbation theory, all fields g are written in the form 
9 = ffo + ffiJ which is the value of the zeroth- 
order state, and gi is the acoustic perturbation which 
by definition has to be much smaller than g^. For the 
velocity, the value of the zeroth-order state is Ug = 0, and 
thus V = v^. The zeroth-order terms solve the governing 
equations in the zeroth-order state and thus drop out of 
the equations. Keeping only first-order terms, we obtain 
the following first-order equations. 

The first-order continuity equation (I3b|) becomes 


KsdtPl = ■ '*^1: 


( 5 ) 


and likewise, the momentum equation (I4cl) becomes 

= V • [ti- pil], ( 6 a) 

where Ti is given by 


n = Vo 


-I- (Vu^)'^ 



(V-uJl. 


( 6 b) 


Equations and ® determine together with a set of 
boundary conditions the time evolution of the first-order 
acoustic fields Pi and v-^. 


When solving for the periodic state at t —oo, it is 
advantageous to formulate the first-order equations in 
the frequency domain. The harmonic first-order fields 
are all written as gi[r,t) = Re {( 7 i'^('r)e“*“*}, where 
is the complex field amplitude in the frequency domain. 
The first-order frequency-domain equations are derived 
from Eqs. and (153 by the substitution —>■ —iw, 

V • = 0, (9) 

V • - pf l] + iwpguf = 0. (10) 

The steady time-averaged streaming flow is obtained 
from the time-averaged second-order frequency-domain 
equations, where { 9 ^ 2 ) denotes time averaging over one 
oscillation period of the periodic second-order field. The 
time-average of products of two harmonic first-order 
fields gf is given by = 5 Re [( 3 “) * 3 “], 

as in Ref. [^, where the asterisk denotes complex con¬ 
jugation. In the periodic state, the fields may consist 
of harmonic terms and a steady term, and thus all full 
time-derivatives average to zero = 0. The time- 

averaged second-order frequency-domain equations are 
derived from Eqs. ® and (Ha)) . 

V-{vi^)+n,{v{^-Vp{^)=0, (11) 

V-[(r“)-(p^'^)l-PoW)]= 0 . ( 12 ) 


D. Second-order time-domain eqnations 


F. Acoustic energy and cavity Q-factor 


Moving on to second-order perturbation theory, we 
write the fields sls g = g^ + gi + g 2 , with g^ and ^2 de¬ 
pending on both time and space. For simplicity and in 
contrast to Ref. m , we do not include perturbations in 77 
and r]^. This will cause the magnitude of the streaming to 
be slightly off, as does the adiabatic approximation, how¬ 
ever the qualitative behavior is not expected to change. 
The second-order time-domain continuity equation (I3bl) 
becomes 

i^s9tP2 =■'^ 2 -KsVi-Vpj_, (7) 


The total acoustic energy of the system in the time 
domain E^{t) and in the frequency domain ^^^^(oo)) is 
given by 





+ 


dF, 


(13a) 




(pfpf) + 


,.fd 


) dF. 
(13b) 


and the momentum equation (I4cl) takes the form 
Po^tV2 = -Pidt^i + V ■ [7-2-7721] - Poi'^i-'^)'^!^ (8a) 


where T 2 is given by 

T 2 = 7?g Vu 2 -f (VU2)'^ 



[ b 2 1 

-f 

Vo - 3 % 


(V • V 2 ) 1. ( 8 b) 


Using Eq. o in the form = p^n^Pi and the first-order 
momentum equation (I 6 al) . we rewrite Eq. (I 8 al) to 


( 8 c) 


Po^t^2 ='^ ■ ["^2 - P2 1 - l^sPl^l + \i^sPI 1] 

-f -Ti -po(Ui • V)Ui. 


This particular form of the second-order momentum 
equation is chosen to minimize numerical errors as de¬ 
scribed in Section UlI Al 


Moreover, the time derivative of is 


= / a 


1 2,1 2 
+ 2^0«1 


dF 


= f [^,p,d,p,+p,v,-dtv,] dF 
Jv 

■ ■Wi ■ (ti-Pi 1)1 -f-dF, (14a) 


where we have used Eqs. dH) and (j^ . Applying Gauss’s 
theorem on the first term in Eq. (|14ap . we arrive at 

dtE^c = / K ■ (U -Pil)j -ndA- [ Vui : dF 
Ja'- J Jv 

= Fpump — This, (14t>) 
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where Ppump is the total power delivered by the forced 
vibration of the sidewalls, and Pdis is the total power 
dissipated due to viscous stress. The quality factor Q of 
a resonant cavity is given by 


Q 


2tt 


Energy stored 
Energy dissipated per cycle 



(15) 


G. Summary of theory 

Throughout this paper we refer to two kinds of so¬ 
lutions of the acoustic energy and velocity fields: un¬ 
steady non-periodic solutions obtained from Eqs. m-m 
and steady periodic solutions obtained from Eqs. ®- 
(ini)- When presenting the unsteady non-periodic solu¬ 
tions, they are often normalized by the steady periodic 
solution, to emphasize how close it has converged towards 
this solution. 


(a) 




Fixed hard wall 


(c) 



h 
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III. NUMERICAL MODEL 

The numerical scheme solves the governing equations 
for the acoustic field inside a water domain enclosed by 
a two-dimensional rectangular microchannel cross sec¬ 
tion. The vibrations in the surrounding chip material and 
piezo transducer are not modeled. The water domain is 
surrounded by immovable hard walls, and the acoustic 
field is excited by oscillating velocity boundary condi¬ 
tions, representing an oscillating nm-sized displacement 
of the walls. A sketch of the numerical model is shown in 
Eig.IlJa). We exploit the symmetry along the horizontal 
center axis z = 0, reducing our computational domain 
by a factor of two. The system is also symmetric about 
the vertical center axis y = 0, however, our attempts to 
use this symmetry introduced numerical errors, and con¬ 
sequently it was not exploited in the numerical model. 
The model used to calculate the steady streaming flow 
in the time-periodic case is a simplification of the model 
presented in Ref. , whereas the model used to solve 
the time-dependent problem is new. 


FIG. 1. (Color online) (a) Sketch of the rectangular compu¬ 
tational domain in the j/z-plane representing the upper half 
of a rectangular cross section of a long straight microchan¬ 
nel of width w = 380 pm and height h = 160 pm as in [^ . 
The thick arrows indicate in-phase oscillating velocity actu¬ 
ation at the left and right boundaries, (b) The three black 
points indicate positions at which the velocity components 
(gray arrows), defined in Eq. (1291) . are probed, (c) Sketch of 
the spatial mesh used for the discretization of the physical 
fields, (d) A zoom-in on the mesh in the upper left corner. 


average of the second-order pressure is enforced by a La¬ 
grange multiplier. Eor the time-domain simulations we 
use the generalized alpha solver [2^, setting the alpha 
parameter to 0.5 and using a fixed time step At. Fur¬ 
thermore, to limit the amount of data stored in Comsol, 
the simulations are run from Matlab and long time¬ 
marching schemes are solved in shorter sections by Com¬ 
sol. Comsol model files and Matlab scripts are provided 
in the Supplemental Material [2^. 


A. Governing equations B. Boundary conditions 


The governing equations are solved using the commer¬ 
cial software Comsol Multiphysics [2^ based on the fi¬ 
nite element method [1^. To achieve greater flexibil¬ 
ity and control, the equations are implemented through 
mathematics-weak-form-PDE modules and not through 
the built-in modules for acoustics and fluid mechanics. 
The governing equations are formulated to avoid eval¬ 
uation of second-order spatial derivatives and of time- 
derivatives of first-order fields in the second-order equa¬ 
tions, as time-derivatives carry larger numerical errors 
compared to the spatial derivatives. To fix the numeri¬ 
cal solution of the second-order equations, a zero spatial 


The acoustic cavity is modeled with stationary hard 
rigid walls, and the acoustic fields are exited on the side 
walls by an oscillating velocity boundary condition with 
oscillation period Iq and angular frequency w. 


The symmetry of the bottom boundary is described by 
zero orthogonal velocity component and zero orthogonal 
gradient of the parallel velocity component. The explicit 
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boundary conditions for the first-order velocity become 


as just the Courant number 


top: 

'<^yl = 0 , 

Vzl = 0 , 

(17a) 

bottom: 

II 

0 

^zl = 0 ^ 

(17b) 

left-right: 

0 

II 

v^i = 0 . 

(17c) 


The boundary conditions on the second-order velocity 
are set by the zero-mass-flux condition n • pu = 0 on 
all boundaries, as well as zero parallel velocity compo¬ 
nent on the top, right and left wall boundaries, and zero 
orthogonal derivative of the parallel component of the 
mass flux on the bottom symmetry boundary. The ex¬ 
plicit boundary conditions for the second-order velocity 
become 


top: 

Vy2 = 0, 

D 

to 

II 

0 

(18a) 

bottom: 

9z{Po^y2+ PlVyl) = 0, 

to 

II 

0 

(18b) 

left-right: 

P0'^y2 + Pl^yl = 0> 

= 0- 

(18c) 


C. Spatial resolution 


The physical fields are discretized using fourth-order 
basis functions for Vi and V 2 and third-order basis func¬ 
tions for Pi and P 2 - The domain shown in Fig. [Ija) is 
covered by basis functions localized in each element of 
the spatial mesh shown in Fig. [TJc). Since the stream¬ 
ing flow is solved in the time domain, the computational 
time quickly becomes very long compared to the compu¬ 
tational time of solving the usual steady streaming flow. 
Thus we have optimized the use of precious few mesh ele¬ 
ments to obtain the best accuracy of the solution. We use 
an inhomogeneous mesh of rectangular elements ranging 
in size from 0.16 pm at the boundaries to 24 pm in the 
bulk of the domain. The convergence of the solution g 
with respect to a reference solution was considered 
through the relative convergence parameter C{g) defined 
in Ref. [ll| by 


C{g) 


\ 


j { 9 - 9reif dydz 

J (ffref)^ dydz 


(19) 


In Ref. [HI, 0 ( 9 ) was required to be below 0.001 for 
the solution to have converged. The solution for the 
steady time-averaged velocity (^^‘^(oo)), calculated with 
the mesh shown in Fig.[Tl(c) and[Tl(d), has C = 0.006 with 
respect to the solution calculated with the fine triangular 
reference mesh in Ref. [ll|, which is acceptable for the 
present study. 


D. Temporal resolution 

The required temporal resolution for time-marching 
schemes is normally determined by the Courant- 
Friedrichs-Lewy (CFL) condition [ 2 ^, also referred to 


CFL = ^ < CFL^ax, (20) 

where At is the temporal discretization and Ar is the 
spatial discretization. This means that the length over 
which a disturbance travels within a time step At should 
be some fraction of the mesh element size, ultimately 
ensuring that disturbances do not travel through a mesh 
element in one time step. A more accurate interpretation 
of the CFL-condition is that it ensures that the error on 
the approximation of the time-derivative is smaller than 
the error on the approximation of the spatial-derivatives. 
Consequently, the value of CFTmax depends on the spe¬ 
cific solver and on the order of the basis functions. For 
fourth-order basis functions and the generalized alpha 
solver, Ref. [1^ reports a value of CFL^^x = 0-05, which 
is an empirical result for a specific model. Due to the in¬ 
homogeneity of the mesh, two values for the upper limit 
for the temporal resolution can be calculated based on 
Eq. (EUl) : At = 8 X 10“^° ns ~ tg/600 for the bulk mesh 
size of 24 pm and At = 5 x 10“^^ ns ~ tg/95000 for the 
boundary mesh size of 160 nm. 

To determine a reasonable trade-off between numeri¬ 
cal accuracy and computational time, we study the con¬ 
vergence of the transient solution towards the steady 
solution for different values of the temporal resolution 
tg/At. The acoustic energy is shown in Fig. [2Ka) 

for different values of At and normalized by the steady 
time-averaged energy (i?^^(oo)) of the frequency-domain 
calculation, and it is thus expected to converge to the 
unity for long times. In Fig.[2Kb), E^^{1000tQ)/ 
is plotted versus the temporal resolution 1^/At, which 
shows how the accuracy of the time-domain solution in¬ 
creases as the temporal resolution is increased. In all 
subsequent simulations we have chosen a time step of 
At = tg/256, the circled point in Fig. [2Kb), for which the 
time-domain energy converge to 99.4% of the energy of 
the steady calculation. The chosen value for the time step 
is larger than the upper estimate tg/600 of the necessary 
At based on the CFL-condition. This might be because 
our spatial domain is smaller than the wavelength, and 
consequently a finer spatial resolution is needed, com¬ 
pared to what is usually expected to spatially resolve a 
wave. 

We have noted that the fastest convergence is obtained 
when actuating the system at its (numerically deter¬ 
mined) resonance frequency /j-gg. When shifting the ac¬ 
tuation frequency half the resonance width 5 A/ away 
from /j.gg, the energy E^^(t) for At = tg/256 converged 
to only 95% of the steady value (i?^g(oo)) (calculated in 
the frequency domain), thus necessitating smaller time 
steps to obtain reasonable convergence. 

The computations where performed on a desktop PC 
with Intel Xeon CPU X5690 3.47 GHz 2 processors, 64- 
bit Windows 7, and 128 GB RAM. The computations 
took approximately one hour for each time interval of 
width lOOtg with At = tg/256, and the computational 
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FIG. 2. (Color online) Numerical convergence and tempo¬ 
ral resolution, (a) Graphs of the build-up of acoustic energy 
Ea,c{t) in the time-domain simulations calculated with differ¬ 
ent fixed time steps At. The energy of the time-domain sim¬ 
ulations is normalized with respect to the energy (^E^{oo)') 
of the steady solution in the frequency domain, and should 
thus converge towards unity. In all simulations the actua¬ 
tion frequency equals the resonance frequency discussed in 
Section IIV Al (b) Acoustic energy i!lac(1000 to) at t = 1000 tpi 
normalized by (i?ac(oo)), and plotted versus the temporal res¬ 
olution tfj/At of the oscillation. The inset is a semilog plot of 
the relative deviation of i?ac(1000 to) from The cir¬ 

cled point in each graph indicates the time step At = to/256 
used in all subsequent simulations. 


time was not limited by RAM, as only less than 2 GB 
RAM was allocated by Comsol for the calculations. 


IV. ONSET OF ACOUSTIC STREAMING 

In this section the fluid is initially quiescent. Then, at 
time t = 0, the oscillatory velocity actuation is turned on, 
such that within the hrst oscillation period its amplitude 
increases smoothly from zero to its maximum value 
which it maintains for the rest of the simulation. We 
study the resulting build-up of the acoustic resonance 
and the acoustic streaming flow. 


A. Resonance and build-up of acoustic energy 

To determine the resonance frequency, the steady 
acoustic energy Eq. (I13b|) was calculated for 

a range of frequencies based on the frequency-domain 
equations (l^- dTUl) . In Fig. |3] the numerical results (cir- 



FIG. 3. (Golor online) Resonance curve and build-up of 
acoustic energy, (a) The numerical acoustic energy density 
(Eac(oo))/V (circles) for different frequencies of the bound¬ 
ary actuation and a Gaussian fit (full line) to the numerical 
data, /rea is the fitted resonance frequency at the center of the 
peak, while /ideal is the frequency corresponding to matching 
a half-wavelength with the channel width. The inset shows 
the numerical build-up of the acoustic energy (full line) for 
actuation at the resonance frequency, uj = 27 r/res, along with 
the analytical prediction Eq. (12311 (dashed line) for a single 
harmonic oscillator with the same resonance frequency and 
quality factor Q — Aes/A/. 


cles) are shown together with a Gaussian fit (full line), 
while the inset exhibits the fitted resonance frequency 
/j.gg, the full width A/ at half maximum, and the quality 
factor Q = 

The build-up of the acoustic energy in the cavity is 
well captured by a simple analytical model of a sin¬ 
gle sinusoidally-driven damped harmonic oscillator with 
time-dependent position x{t), 

H ^ T H T* 1 

— -b 2ra;o— -b ojIx = —Fq sin(a;t). (21) 
aF at m 

Here, T is the non-dimensional loss factor, uJq is the res¬ 
onance frequency of the oscillator, is the amplitude 
of the driving force divided by the oscillator mass, and 
oj is the frequency of the forcing. The loss factor is re¬ 
lated to the quality factor by T = 1/{2Q), and in the 
underdamped case T < 1, the solution becomes 


x{t) = A 


sin(wt -b 4>) 


uje ^‘^0* . / r -- 

- A -F2 sinfVl -T^WotT 

WqVI - V 


( 22 ) 


The amplitude A and the phase shift ip between the forc¬ 
ing and the response are known functions of ■^, Wg, w, 
and r, which are not relevant for the present study. From 
Eq. (|22ll we obtain the velocity dx/dt, leading to the total 
energy E of the oscillator. 


E = -b \m 



(23) 
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Based on Eqs. (l22l) and (l23t . the characteristic timescale 
Tg for the build-up of the acoustic energy is found to be 


2ra;n 


R. 

Wo' 


(24) 


the perturbation expansion becomes invalid for smaller 
values of than previously expected. 

In the transient regime we cannot Fourier decompose 
the velocity field. Instead, we propose a decomposition 
using envelope functions inspired by Eq. (EH), 


The build-up of the energy in the single harmonic os¬ 
cillator, calculated at w = Wq with E = 1.20 x 10“^, 
is shown in the inset of Fig. |3] together with the build¬ 
up of acoustic energy of the microfluidic channel 

solved numerically at resonance, w = 27r/,,gg. The ana¬ 
lytical and numerical results are in good agreement, and 
we conclude that the build-up of acoustic energy in the 
channel cavity can be modeled as a single harmonic os¬ 
cillator. The energy builds up to 95% of its steady value 
in about 500 tg ~ 8 t^. 

B. Decomposition of the velocity field 


u(r, t) = Vi (r, t) sin(wt) -I- t) sin(2a;t) -I- V 2 {r, t). 

(29) 

Here, the amplitudes are slowly varying in time compared 
to the fast oscillation period tg. We can no longer sep¬ 
arate and ^2 before solving the second-order time- 
dependent equations ED and dl])- To obtain the time- 
dependent magnitude of the quasi-steady streaming ve¬ 
locity mode ^2, we need to choose a good velocity probe, 
and we thus form the unsteady time-average of V 2 {r,t)^ 

(« 2 (»’: 0 )=/ V2{r,t')dt'. ( 30 ) 

Jt-tg/2 


The task of calculating the build-up of the acoustic 
streaming flow is a multi-scale problem, because the am¬ 
plitude of the oscillating acoustic velocity field is several 
orders of magnitude larger than the magnitude of the 
streaming flow. This is indeed the very reason that we 
can apply the perturbation expansion 


The time-averaging is done with a fifth-order Romberg 
integration scheme using data points with a uniform 
spacing of Iq/IG in the time interval of width tg. 

C. Steady and unsteady streaming flow 


V = V^ +V2, 


(25) 


and decompose the non-linear governing equations into 
a set of linear first-order equations and a set of second- 
order equations. However, there is also another level of 
difference in velocity scaling. In the purely periodic state, 
the velocity can be Fourier decomposed as 

v{r,t) = Vi{r) sin(wf) -I- ^^“(r) sin(2wf) -I- ^^(t'), (26) 


where u“(r) is the steady amplitude of the first-order 
harmonic component, V 2 ^{r) is the steady amplitude 
of the second-order frequency-doubled component, and 
V 2 {r) is the magnitude of the second-order steady veloc¬ 
ity component referred to as the acoustic streaming flow. 
The orders of magnitude of the three velocity components 
in the periodic state are given by 


Qv. 


be’ 


1 ?^ 

^2 


QS 


be 


^2° 


Q'^v 


be 


(27) 


The order of Vi is derived in the one-dimensional acoustic 
cavity example presented in Ref. [s^, the order of V 2 
is given by the well-known Rayleigh theory, while the 
order of is a new result derived in Appendix A. The 
magnitude of is a factor of Q larger than what is 
expected from dimensional analysis of the second-order 
equation (1^ . Consequently, the criterion |u2| |u;^| for 

the perturbation expansion becomes 




(28) 


which is more restrictive than the usual criterion based on 
the first-order perturbation expansion, <C c^. Thus, 


In this section we compare the unsteady time- 
averaged second-order velocity field (^2(7',^)), from 
the time-domain simulations, with the steady time- 
averaged second-order velocity field (^^‘^(r,00)), from 
the frequency-domain simulation. Figure SDs-) (b) 
each shows a snapshot in time of the transient and 
V 2 , respectively. For V 2 {r,t), the oscillatory component 
u|‘^(r, t) sin(2ajt) dominates, as it is two orders of mag¬ 
nitude larger than the quasi-steady component V 2 (r,t). 
However, at late times, here t = 3000 tg, the amplitude 
has converged, and in (v 2 {r,t)'^ the oscillatory 
component average to zero and only the quasi-steady 
component remains. 

The unsteady time average (v 2 {r,t)) evaluated at 
t = 3000 tg is shown in Fig.|4Kc), exhibiting a single flow 
roll, in agreement with the classical Rayleigh stream¬ 
ing flow. In Fig. EDd) is shown the steady (^^‘^(oo)) 
from the frequency-domain simulation. Figure SDc) and 
EKd) use the same color scaling for the velocity magni¬ 
tude, to evaluate the convergence of the unsteady stream¬ 
ing flow (1)2(3000 tg)) towards the steady streaming flow 
(uf^(oo)), and the two solutions agree well both quali¬ 
tatively and quantitatively. The convergence parameter 
C, Eq. (fTOl) . of (7)2(3000tg)) with respect to (^^‘^(oo)) 
is C = 0.01, and if we multiply (7)2) by a free factor, 
taking into account that (7)2) has not fully converged at 
t = 3000 tg, the convergence parameter can be reduced 
to C = 0.008. The remaining small difference between 
the unsteady (7)2(3000 tg)) and the steady (00)) is 
attributed to the finite temporal resolution of the time 
marching scheme. We can thus conclude that the time- 
domain streaming simulation converges well towards the 





frequency-domain simulation, and this constitutes the 
primary validation of the unsteady non-periodic simu¬ 
lations. 


(a) I 







(c) 


0 

h 

2 



0 

(d)| 





0 





FIG. 4. (Color online) (a) Snapshot of the oscillatory first- 
order velocity field Vi (vectors) and its magnitude [color plot 
ranging from 0 m/s (black) to 0.7 m/s (white)] at t = 3000 Iq. 
(b) Snapshot of the oscillatory second-order velocity field 
(vectors) and its magnitude [color plot ranging from 0 m/s 
(black) to 0.02 m/s (white)[ at t = 3000 tg. (c) Snapshot of the 
unsteady time-averaged second-order velocity field (U 2 ) (vec¬ 
tors), Eq. (I30II . and its magnitude [color plot ranging from 0 
mm/s (black) to 0.1 mm/s (white)] at t = 3000 Iq. (d) Steady 
time-averaged second-order velocity field (u|‘^(oo)) (vectors), 
Eqs. (Illll and inj, and its magnitude [color scaling as in (c)]. 
In both the time-domain and the frequency-domain simula¬ 
tions the parameters of the oscillating velocity boundary con¬ 
dition was LJ = 27r/res and = tod, with wall displacement 
d = 1 nm. 


D. Build-up of the velocity field 


To visualize the build-up of the acoustic fields over 
short and long timescales, we have chosen the three point 
probes shown in Fig.[T][b). The oscillating first-order ve¬ 
locity field is probed in the center of the channel (0,0), 
far from the walls in order to measure the bulk ampli¬ 
tude of the acoustic held. The horizontal component of 
the second-order velocity Vy 2 is probed on the horizontal 
symmetry axis at (|;rc, 0), where the oscillatory compo¬ 
nent 1)1“^ has it maximum amplitude. The vertical com¬ 
ponent of the second-order velocity v ^2 probed on the 
vertical symmetry axis at (0, -jh) where the oscillatory 
component is small and of the same order as the 
quasi-steady component v^., making the unsteady time- 
averaged second-order velocity at this point a good probe 
for the quasi-steady streaming velocity. 

In Fig. [5] is shown the build-up of the velocity probes 
(a-c) and their time-averages (d-f) for the hrst 20 os¬ 
cillations. The thick lines are the oscillating velocities 
while, the thin lines are the envelopes of the oscillations. 
Already within the first 20 oscillation periods we see in 
Fig. [SKf) the build-up of a quasi-steady velocity com¬ 
ponent. The unsteady time-averaged horizontal veloc¬ 
ity (v)’ Fig-lHe), is still primarily oscillatory, showing 
that for this probe the oscillatory component is much 
larger than the quasi-steady component t>2. 

The temporal evolution of the velocity probes on the 
longer time scale up to t = 1500 Iq is shown in Fig. [51 
In Fig. m a) and (b) the amplitudes of the oscillatory 
first- and second-order velocity components are seen to 
stabilize around t = 700 tp ^ 10The steady ampli¬ 
tudes of the velocity probes Fig. [S] agree with the the¬ 
oretical predictions of Eq. (123, yielding orders of mag¬ 
nitude Vi/Cg ^ 3 X 10“"^ (Fig. [6l[a)), jc^ ~ 5 x 10“^ 
(Fig. HDb)), and v^jc^ ~ 1 x 10“’’ (Fig. [^e) andlHKf))- 
The time-average of Vyi tends to zero for long times 
as it is purely oscillatory, whereas the time-average of 
Vy 2 tends to the magnitude of the quasi-steady compo¬ 
nent V 2 , because the large but now steady oscillatory 
component average to zero. The dashed lines in 
Fig. [He) and (f) represent the magnitude of the steady 
time-averaged second-order velocity (^^‘^(oo)) from the 
frequency-domain simulation. 


V. ACOUSTIC STREAMING GENERATED BY 
PULSED ACTUATION 


In the following we study the effects of switching the 
oscillatory boundary actuation on and off on a timescale 
much longer than the oscillation period tg in either 
single- or multi-pulse mode. The aim is to investigate 
whether such an approach can suppress the influence of 
the streaming flow on suspended particles relative to that 
of the radiation force. 
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FIG. 5. (Color online) Velocity probes for the initial time interval 0 < t < 20to- (S'-c) probes for the first- and second-order 
velocity (a) nyi(0,0), (b) Vy 2 {w/A,Q), and (c) Vz2(Q,h/4). (d-f) Running time-average Eq. (1301) on an interval one oscillation 
period wide of the velocity probes in (a-c). The thick lines show the oscillating velocity probes while the thin lines emphasize 
the envelopes of the oscillations. 



FIG. 6. (Color online) The velocity probes from Fig. [5] but now extended to the long time interval 0 < t < 1500 to, showing 
the convergence towards a periodic state. The dashed lines in (e) and (f) indicate the magnitude of the steady time-averaged 
second-order velocity from the frequency-domain simulation Eqs. m and HU. 
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TABLE II. Characteristic timescales. The values are obtained 
by using the kinematic viscosity v = rj/Po = 8-93 x 10“^ m^/s 
(Table HI), the Q-factor Q = 416 (Fig. [Sj, and the channel 
height h — 160 jim (Fig. [TJ . 


Timescale 

Expression 

Value 


Oscillation time 

^0 

5.1 X lO"’’ s R 

8 1 to 

Resonance 




relaxation time 

-^E = ^to 

3.4 X 10 ® s R 

s 66 to 

Momentum 




diffusion time 

2i^ ( 8 ) 

2.8 X 10 ^ s R 

s 558 to 


A. Single-pulse scaling analysis 

A striking feature of Fig. El is the separation of 
timescales between the roughly exponential build-up of 
the acoustic resonance in Fig. E^a) and of the stream¬ 
ing flow in Fig. EIO- It appears that the resonance, and 
hence the acoustic radiation force on a suspended par¬ 
ticle, is fully established almost ten times faster than 
the streaming flow and the resulting drag force on a sus¬ 
pended particle. To investigate this further, we look at 
the scaling provided by the three timescales relevant for 
the problem of transient acoustic streaming, all listed 
in Table m the oscillation time of the acoustic wave, 
the resonance relaxation time of the acoustic cavity, 
and the momentum diffusion time governing the quasi¬ 
steady streaming flow. 

The momentum diffusion time is , where 

1/ = ^ is the kinematic viscosity, and is approxi¬ 
mately half the distance between the top boundary and 
the center of the streaming flow roll. Inserting the rele¬ 
vant numbers, see TablelU we indeed find that te ~ 66 Iq 
is much faster than Ri 558tg. However, this separa¬ 
tion in timescales does not guarantee a suppression of 
streaming relative to the radiation force. One problem is 
that the streaming is driven by the shear stresses in the 
boundary layer, and these stresses builds up much faster 
given the small thickness of the boundary layer. This we 
investigate further in the following subsection. Another 
problem is that the large momentum diffusion time 
implies a very slow decay of the streaming flow, once it is 
established. The latter effect, we study using the follow¬ 
ing analytical model. Consider a quantity / (streaming 
velocity or acoustic energy), with a relaxation time r and 
driven by a pulsed source term P of pulse width 
rate of change of / is equivalent to Eq. (I14bl) . 

dJ = P--f, 

r 

i/o, for 0 < t < tp^, 

0, otherwise, 

where ^/g is a constant input power. This simplified 
analytical model captures the roughly exponential build¬ 
up and decay characteristics of our full numerical model. 



tpw The 

(31a) 

(31b) 


and allows for analytical studies of the time-integral of 
f{t). For a final time t > we Hnd 


[ f{t’) dt' = /gtp^ - foT 

Jo 


3-v(*-V) - 


(32) 


From this we see that when t ^ t +the time-integral 
of f(t) is approximately /otpw and not dependent on the 
relaxation time r. Consequently, if both the acoustic en¬ 
ergy and the acoustic streaming can be described by ex¬ 
ponential behavior with the respective relaxation times 
Te and the ratio of their time-integrated effects is 
the same whether the system is driven by a constant ac¬ 
tuation towards their steady time-periodic state or by 
a pulsed actuation with pulse width tp.^^. This simpli- 
hed analytical model indicates that there is little hope 
of decreasing acoustic streaming relative to the acoustic 
radiation force by applying pulsed actuation, in spite of 
the order of magnitude difference between the relaxation 
times for the acoustic energy and the streaming. 


B. Single-pulse numerical analysis 


We investigate the features of pulsed actuation more 
detailed in the following by numerical analysis. In Fig. [7] 
is shown the temporal evolution of the total acoustic en¬ 
ergy Ih® magnitude of the acoustic streaming 

flow (vgtr) for the three cases: (i) the build-up towards 
the periodic state, (ii) a single long actuation pulse, and 
(iii) a single short actuation pulse. The magnitude of the 
acoustic streaming is measured by the unsteady time- 
averaged velocity probe 

(Vstr) = {Vz2id,^h)), (33) 

and the unsteady energy and streaming probes obtained 
from the time-domain simulation are normalized by 
their corresponding steady time-averaged values from the 
frequency-domain simulation. 

We introduce the streaming ratio x to measure the in¬ 
fluence of streaming-induced drag on suspended particles 
relative to the influence of the acoustic radiation force 
for the unsteady time-domain solution, in comparison to 
the periodic frequency-domain solution. To calculate the 
relative displacement As of particles due to each of the 
two forces, respectively, we compare their time integrals. 
Since the radiation force scales with the acoustic energy 
density, we define the streaming ratio x(t) as 


(^str(^')) , 

^ Jo (^str(oo)) _ Asg, 

f )) , As^.,^g 

Jo (E«(oo)) As“g 


(34) 


where Asgtr and Asrad are the total particle displace¬ 
ments in the time from 0 to t due to the streaming- 
induced drag force and the acoustic radiation force, re¬ 
spectively. In the periodic state y = 1, and to obtain 
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FIG. 7. (Color online) Acoustic energy (l?ac(l))/(^'ic(oo)) 
Eq. (IT 3 I 1 (light green), streaming velocity (wstr(l))/(^^str(oo)) 
Eq. (1331) (medium purple), and streaming ratio x(l) Eq. (1341) 
(dark brown, right ordinate axis). The gray background indi¬ 
cates the time intervals where the actuation is turned off. 
(a) Constant actuation for 0 < t < 3000to- (b) Actua¬ 
tion on for 0 < t < 200 to followed by no actuation for 
200 to < t < 1000 to. (c) Actuation on for 0 < t < 30to 
followed by no actuation for 30to < t < 1000 to. 


radiation force-dominated motion of smaller particles, 
we need to achieve a smaller value of x- Obtaining a 
value of X = 0.8 at time tgndJ implies that the ratio of 
the relative displacement due to the streaming-induced 
drag force and the radiation force for the time interval 
0 < t < is 20% lower than in the periodic state, cor¬ 
responding to a 20% reduction of the critical particle size 
for acoustophoretic focusing, defined in Ref. 0, assum¬ 
ing the particles can be focused during the time interval 
0 < i < iend- 

Figure[7Ka) shows ('*^str) S'lid x during the build¬ 

up towards the periodic state, x approaches unity slower 
than ('Cgtr) because x is an integration of the stream¬ 
ing and radiation contributions, whereas probes the 
instantaneous magnitude of the streaming flow. Figure 



t/tg 


FIG. 8. (Color online) The same probes as in Fig. [7] but for 
the following pulsed actuation schemes: (a) actuation is on 
for 500 to followed by no actuation for 500 to repeatedly, (b) 
actuation is on for 200 to followed by no actuation for 200 to 
repeatedly, and (c) actuation is on for 30 to followed by no 
actuation for 210 to repeatedly. 


[T^b) and[7Kc) show (E^^), (fstr): x when the actua¬ 

tion is turned off at t = 200 Iq and t = 30 tp, respectively. 
When the actuation is turned off, (ifac) decays faster 
than (ugjj.) and thus x begins to increase more rapidly, 
reaching x = 0.8 around t = 1000 tg in both cases. From 
the results shown in Fig. [7] it does not seem advantages to 
turn off the actuation, as this only causes x to increase 
faster than for constant actuation. Figure [7Kc) further 
shows that when the actuation is turned off, {E^') im¬ 
mediately begins to decay, whereas (r’gtr) continues to 
increase for some time, due to the present acoustic en¬ 
ergy in the system that still provides a driving force for 
the streaming flow. 
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C. Multi-pulse numerical analysis 

From the single pulse results shown in Fig.[7]there is no 
indication of any optimum for the pulse duration or rep¬ 
etition period, and in general it provides little hope that 
pulsed actuation should lead to lower values of %. Fig¬ 
ure [8] shows (i?ac)i (^str)i X for three pulsed schemes 
with pulse duration 500 tp, 200 tg, and 30 tg and pause du¬ 
ration 500 tg, 200 tg, and 210 tg, respectively. For all three 
pulsed schemes, x increases faster than for the constant 
actuation Fig. [TJa), thus not indicating any increased 
suppression of the streaming. 


VI. DISCUSSION 

Solving numerically the time-dependent problem of the 
acoustic cavity and the build-up of acoustic streaming, 
presents new challenges, which are not present in the 
purely periodic problem. Firstly, the numerical conver¬ 
gence analysis now involves both the spatial and tempo¬ 
ral resolutions. This we addressed in a sequential pro¬ 
cess by first analyzing the spatial mesh with the periodic 
frequency-domain solution, and thereafter doing a thor¬ 
ough convergence analysis with respect to the temporal 
resolution. Secondly, the convergence of the transient 
solution towards the periodic state was poor for actu¬ 
ation frequencies away from the resonance frequency of 
the system. This makes off-resonance simulation compu¬ 
tationally costly, as it requires a better temporal reso¬ 
lution, and it complicates comparison of simulations at 
resonance with simulations off resonance. Thirdly, small 
numerical errors accumulate during the hundred thou¬ 
sand time steps taken during a simulation from a quies¬ 
cent state to a purely periodic state. These errors need 
to be suppressed by the numerical time-domain solver, 
which in the generalized-alpha solver is done through the 
alpha parameter. Simulation with higher temporal res¬ 
olution required lower values of the alpha parameter to 
have more suppression of accumulated numerical errors. 

The model system used in this study is a simplifica¬ 
tion of an actual device. The vibration of only the side 
walls, and not the top and bottom walls, stands in con¬ 
trast to the physical system, in which the whole device 
is vibrating in a non-trivial way, difficult to predict, and 
only the overall amplitude and the frequency of the ac¬ 
tuation is controlled experimentally. Furthermore, our 
model only treats the two-dimensional cross section of a 
long straight channel, whereas experimental studies have 
shown that there are dynamics along the length of the 
channel [l^. Nevertheless, successful comparison, both 
qualitatively and quantitatively, have been reported be¬ 
tween the prediction of this simplified numerical model 
and experimental measurements of Rayleigh streaming 
in the cross sectional plane of a microchannel 0 , which 
makes it reasonable to assume that the time-dependent 
simulations also provide reliable predictions. 

It is also important to stress that our model only de¬ 


scribes the fluid and not the motion of the suspended 
particles. Integrating the forces acting on the particles 
becomes vastly more demanding when the streaming flow 
is unsteady, because the drag forces from the oscillating 
velocity components and do average out, as they 
do in the case of a purely time-periodic state. To in¬ 
clude this contribution in the particle tracking scheme, 
the forces on the particles need to be integrated with a 
time step of a fraction of the oscillation period, which 
makes the solution of particle trajectories over several 
seconds a very demanding task using brute-force integra¬ 
tion of the equations of motion. 

Our analysis of the pulsed actuation schemes showed 
that the slow decay of the streaming flow makes pulsa¬ 
tion inefficient in reducing the streaming-induced drag 
force compared to the radiation force. Such a reduction 
may, however, be obtained by a rapid switching between 
different resonances each resulting in similar radiation 
forces but different spatial streaming patterns which on 
averages cancel each other out, thus fighting streaming 
with streaming. An idea along these lines was presented 
by Ohlin et al. Ref. (^ . who used frequency sweeping 
to diminish the streaming flows in liquid-filled wells in a 
multi-well plate for cell analysis. However, the prediction 
of particle trajectories under such multi-resonance con¬ 
ditions requires an extensive study as described above. 

Experimentally, the use of pulsed actuation to decrease 
streaming flow has been reported by Hoyos et al. Ref. 

However, this study is not directly comparable to 
our analysis, as we treat the build-up of Rayleigh stream¬ 
ing perpendicular to the pressure nodal plane, whereas 
Hoyos et al. study the streaming flow in this plane. Such 
in-nodal-plane streaming flows have been studied numer¬ 
ically by Lei et al. , though only with steady actu¬ 

ation. The contradicting results of our theoretical study 
and the experimental study of Hoyos et al. may thus rely 
on the differences of the phenomena studied. 


VII. CONCLUSION 

In this work, we have presented a model for the 
transient acoustic fields and the unsteady time-averaged 
second-order velocity field in the transverse cross- 
sectional plane of a long straight microchannel. The 
model is based on the usual perturbation approach for 
low acoustic field amplitudes, and we have solved both 
first- and second-order equations in the time domain for 
the unsteady transient case as well as in the frequency 
domain for the purely periodic case. This enabled us to 
characterize the build-up of the oscillating acoustic fields 
and the unsteady streaming flow. 

Our analysis showed that the build-up of acoustic en¬ 
ergy in the channel follows the analytical prediction ob¬ 
tained for a single damped harmonic oscillator with si¬ 
nusoidal forcing, and that a quasi-steady velocity compo¬ 
nent is established already within the first few oscillations 
and increases in magnitude as the acoustic energy builds 
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up. We have also found that for a resonance with quality 
factor Q, the amplitude of the oscillatory second-order 
velocity component is a factor of Q larger than what is 
expected from dimensional analysis, which results in a 
more restrictive criterion for the validity of the perturba¬ 
tion expansion, compared to the usual one based on the 
first-order perturbation expansion. 

Furthermore, contrary to a simple scaling analysis of 
the time scales involved in the fast build-up of radia¬ 
tion forces and slow build-up of drag-induced streaming 
forces, we have found that pulsating oscillatory boundary 
actuation does not reduce the time-integrated streaming- 
induced drag force relative to the time-integrated radia¬ 
tion force. As a result, pulsating actuation does not pre¬ 
vent streaming flows perpendicular to the pressure nodal 
plane from destroying the ability to focus small particles 
by acoustophoresis. 
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Appendix A: Amplitude of the second-order 
oscillatory velocity field 

Extending to second order the one-dimensional exam¬ 
ple given in Ref. [s^ , we derive in this appendix the order 
of magnitude of the second-order oscillatory component 
which was stated in Eq. (I2Z1). 

Like (52) denotes time-averaging over one oscillation 
period, Eq. (1501) . and in the periodic state equals the 
zero-order temporal Fourier component of the field, then 
g 2 ^{r) denotes the complex amplitude of the oscilla¬ 
tory second-order mode and is given by the second-order 
Fourier component 

9 r{r) = - dt'. (Al) 

J Jt-T/2 

By using the general formula for the real part 
of any complex number Z, Re[Z] = ^(Z + Z*), 
the product A{r,t)B{r,t) of two oscillating fields 
A{r,t) = Re [Ae““^‘] and B{r,t) = Re [i?e“'“*] can be 
decomposed into a steady component and an oscillatory 


component 


A(t)Bit) = ifAe-" 


from which we introduce the following notation 


(AB) = i Re 


A*B 


(A3) 


{ABf‘^ = \AB, (A4) 

where A and B could be any first-order helds. 

The governing equations for the oscillatory second- 
order component can be derived from Eqs. © 

and (fSl) and in the one-dimensional problem treated in 
Ref. 1^, where the top and bottom walls are not taken 
into account, they become 


-i2uJK^pl‘^ = -dyvf" - (A5a) 

-\ 2 u:py 2 ^ = -dypl^ + (fr? + rf^) 

- (Pi(-iwui))^‘^ - p^iy^dyV^f'^. (A5b) 

Applying the 2w-rule of Eq. (IA4I) and the mass continuity 
Eq. ([5]), the two last terms of Eq. (|A5b(l cancel. Inserting 
Eq. (lASal) into Eq. (lASbl) . the governing equation for 
becomes 

Aklvf^ + (1 - 'AT)d^vl^ + ]^K^dy{v^dyP^) = 0, (A6) 


where F is the non-dimensional bulk damping coeffi¬ 
cient given by F = 2^(| + is the 

wavenumber. For the fundamental half-wave resonance, 
the spatial dependence of the source term dy{vidyPi) is 
sin(2A:Q?/), and the guess for the inhomogeneous solution 
to Eq. (IA6I) thus becomes 


U2 ’ =Csin(2/fco?/). (A7) 

Inserting the inhomogeneous solution Eq. (IAtI) into the 
governing equation (IA6I) . we note that the first term can¬ 
cels with the “1” in the parentheses of the second term, 
and the order of magnitude of the inhomogeneous solu¬ 
tion thus becomes 


,,2cj I 


= c^ 




F3 c. 


Q 




(A8) 


which is the result stated in Eq. (EZl). 
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